
** Replication code for Taylor C. Boas and Amy Erica Smith, "Looks Like Me, Thinks Like Me: Descriptive Representation and Opinion Congruence 
** in Brazil." Latin American Research Review 54, 3 (2019).

** Analysis conducted in StataSE Version 14 on Windows 10 Enterprise.


** NOTE: This file replicates main text Figure 1 and Appendix Figures 2 and 3. 
** It uses data from the Churches North and South project, conducted Fall 2014 by Amy Erica Smith.
** Both the clergy and congregant datasets may be downloaded from: https://dataverse.harvard.edu/dataverse/religion_democracy_brazil


** Set working directory as appropriate

**************************************************************************************************************************************
********************************************* 1. CLERGY SURVEY ANALYSIS **************************************************************

use "ChurchesNorthAndSouth_Brazil2014_Clergy", clear
g juiz = city == "JF"
recode p1_denomination (1=1) (2 3 5 = 2) (4 = 3) (6 = .), g(denom)
	lab def denom 1 "Catholic" 2 "Evangelical/Mainline" 3 "Pentecostal"
	lab val denom denom

********** Figure 3: Frequency of discussion of selected topics among Juiz de Fora clergy
preserve
drop if juiz != 1
recode denom (3=2)
collapse (mean) p26_discuss_racism p24_discuss_abortion p23_discuss_homosexuality p25_discuss_environment ///
	(semean) se_ra=p26_discuss_racism se_ab=p24_discuss_abortion se_gay=p23_discuss_homosexuality se_env=p25_discuss_environment, by(denom)
drop if denom == .
set obs 12
egen denom2 = seq(), f(1) t(3) block(1)
	lab val denom2 denom
	drop denom
	rename denom2 denom
foreach i in p26_discuss_racism p24_discuss_abortion p23_discuss_homosexuality p25_discuss_environment se_ra se_ab se_gay se_env {
	egen `i'mn = mean(`i'),by(denom)
	replace `i' = `i'mn
	drop `i'mn
}
g xaxis = _n
egen dv = seq(), f(1) t(4) block(3)
	lab def dv 1 "Combat Racism" 2 "Abortion" 3 "Homosexuality" 4 "Environment", modify
	lab val dv dv
g dv_val = p26_discuss_racism if dv == 1
	replace dv_val = p24_discuss_abortion if dv == 2
	replace dv_val = p23_discuss_homosexuality if dv == 3
	replace dv_val = p25_discuss_environment if dv == 4
g se_val = se_ra if dv == 1
	replace se_val = se_ab if dv == 2
	replace se_val = se_gay if dv == 3
	replace se_val = se_env if dv == 4
g lb = dv_val-2*se_val
g ub = dv_val+2*se_val
graph twoway (bar dv_val xaxis if denom == 1, lcolor(black) fcolor(gs0)) || ///
			(bar dv_val xaxis if denom == 2, lcolor(black) fcolor(gs10)) || ///
			(rspike ub lb xaxis, lcolor(black)), ///
		legend(order(1 "Catholic" 2 "Evangelical") cols(2) symxsize(*.25) size(medsmall)) ///
		xlabel(1.5 "Racism" 4.5 "Abortion" 7.5 "Homosexuality" 10.5 "Environment", notick) xtitle("") ///
		ylabel(1(1)5) ytitle("Average Reported Frequency (1-5 Scale)," "with 95% Confidence Intervals") ///
		graphregion(fcolor(white) lcolor(white)) saving(clergy_figure.pdf) 
restore

********** Alternate version of Figure 3: all clergy (not only Juiz de Fora)
preserve
recode denom (3=2)
collapse (mean) p26_discuss_racism p24_discuss_abortion p23_discuss_homosexuality p25_discuss_environment ///
	(semean) se_ra=p26_discuss_racism se_ab=p24_discuss_abortion se_gay=p23_discuss_homosexuality se_env=p25_discuss_environment, by(denom)
drop if denom == .
set obs 12
egen denom2 = seq(), f(1) t(3) block(1)
	lab val denom2 denom
	drop denom
	rename denom2 denom
foreach i in p26_discuss_racism p24_discuss_abortion p23_discuss_homosexuality p25_discuss_environment se_ra se_ab se_gay se_env {
	egen `i'mn = mean(`i'),by(denom)
	replace `i' = `i'mn
	drop `i'mn
}
g xaxis = _n
egen dv = seq(), f(1) t(4) block(3)
	lab def dv 1 "Combat Racism" 2 "Abortion" 3 "Homosexuality" 4 "Environment", modify
	lab val dv dv
g dv_val = p26_discuss_racism if dv == 1
	replace dv_val = p24_discuss_abortion if dv == 2
	replace dv_val = p23_discuss_homosexuality if dv == 3
	replace dv_val = p25_discuss_environment if dv == 4
g se_val = se_ra if dv == 1
	replace se_val = se_ab if dv == 2
	replace se_val = se_gay if dv == 3
	replace se_val = se_env if dv == 4
g lb = dv_val-2*se_val
g ub = dv_val+2*se_val
graph twoway (bar dv_val xaxis if denom == 1, lcolor(black) fcolor(gs0)) || ///
			(bar dv_val xaxis if denom == 2, lcolor(black) fcolor(gs10)) || ///
			(rspike ub lb xaxis, lcolor(black)), ///
		legend(order(1 "Catholic" 2 "Evangelical") cols(2) symxsize(*.25) size(medsmall)) ///
		xlabel(1.5 "Racism" 4.5 "Abortion" 7.5 "Homosexuality" 10.5 "Environment", notick) xtitle("") ///
		ylabel(1(1)5) ytitle("Average Reported Frequency (1-5 Scale)," "with 95% Confidence Intervals") ///
		graphregion(fcolor(white) lcolor(white)) 
restore


**************************************************************************************************************************************
********************************************* 2. CONGREGANT SURVEY ANALYSIS **********************************************************

use "ChurchesNorthAndSouth_Brazil2014", clear
clonevar site = igreja if igreja != 7
	replace site = 10+communitylocation if site == .
egen site2 = group(site) 
recode site2 (2 3 4 8 =1) (1 5 6 7 = 2) (9 10 11 12 13 = 3), g(relig)
	lab def relig 1 "Catholic" 2 "Evangelical/Pentecostal" 3 "Non-Church Site"
	lab val relig relig

********** Figure 4: Variance in support for policy issues, by type of community site
global coefficients ""
foreach i in 21 18 19 20 {
tempfile cath evan comm
reg p`i'_ i.site2
	predict p`i'resid, resid
	replace p`i'resid = (p`i'resid)^2
	reg p`i'resid b2.relig 
estimates store estp`i'
	lincomest _cons
	parmest, saving(`evan'_`i', replace) flist(coefficients)
estimates restore estp`i'
	lincomest _cons + 1.relig
	parmest, saving(`cath'_`i', replace) flist(coefficients)
estimates restore estp`i'
	lincomest _cons + 3.relig
	parmest, saving(`comm'_`i', replace) flist(coefficients)
drop p`i'resid _est_estp`i'
}
preserve
dsconcat $coefficients
egen site = seq(), f(1) t(3)
	lab def site 1 "Evangelical Churches" 2 "Catholic Churches" 3 "Community Sites"
	lab val site site
egen dv = seq(), f(1) t(4) block(3)
	lab def dv 1 "Environment" 2 "Gay Marriage" 3 "Abortion" 4 "Combat Racism" 
	lab val dv dv
g yaxis = _n + dv - 1
graph twoway (rspike min95 max95 yaxis if site == 1, horizontal lcolor(black)) ///
			(rspike min95 max95 yaxis if site == 2, horizontal lcolor(gs5)) ///
			(rspike min95 max95 yaxis if site == 3, horizontal lcolor(gs10)) ///
			(scatter yaxis estimate if site == 1, msymbol(O) mcolor(black)) ///
			(scatter yaxis estimate if site == 2, msymbol(T) mcolor(gs5)) ///
			(scatter yaxis estimate if site == 3, msymbol(X) mcolor(black)), ///
	xscale(r(0 )) xlabel(0(.5)3, nogrid) xline(0, lcolor(gs10)) ///
	yscale(r(.25 16)) ylabel(2 "Environment" 6 "Gay Marriage" 10 "Abortion" 14 "Racism", nogrid angle(horizontal) notick) ///
	yline(4 8 12 16, lcolor(gs13) lpattern(dash)) ytitle("") ///
	legend(order(6 "Non-Church" 5 "Catholic" "Churches" 4 "Evangelical" "Churches" ) size(small) col(1) ///
			ring(0) position(1) region(margin(small))) ///
	graphregion(fcolor(white) lcolor(white)) plotregion(lcolor(black)) ///
	xtitle("Mean Squared Residuals, with 95% Confidence Intervals", color(black)) saving(congregation_variance_figure.pdf)
restore


********** Appendix Table with coefficients from Figure 4
preserve
global coefficients ""
foreach i in 21 18 19 20 {
reg p`i'_ i.igreja
	estimates store estp`i'
predict p`i'resid, resid 
	replace p`i'resid = (p`i'resid)^2
reg p`i'resid b2.relig
	estimates store estp`i'resid
}
esttab estp21 estp18 estp19 estp20 using vfr_firststage_reg.tex, ///
		b(3) se(3) nostar replace label nobaselevels mtitles("Environment" "Gay Marriage" "Abortion" "Racism") ///
		addnotes("Coefficients represent results from models regressing policy attitudes on a" ///
		"categorical variable for site. Evangelical Church 1 is the excluded category." ///
		"Standard errors in parentheses.") nonote // table of first stage regressions
esttab estp21resid estp18resid estp19resid estp20resid using vfr_secondstage_reg.tex, ///
		b(3) se(3) nostar replace label nobaselevels mtitles("Environment" "Gay Marriage" "Abortion" "Racism")  ///
		addnotes("Coefficients represent results from models regressing squared residuals from first" ///
		"stage on a categorical variable for type of site. Evangelical church is the excluded" ///
		"category. Standard errors in parentheses.") nonote // table of second stage regressions
restore
	
********** Appendix: analysis reestimated using gamma regression, following Western and Bloome (2009)
foreach i in 21 18 19 20 {
global coefficients ""
tempfile cath evan comm
reg p`i'_ i.site2
	predict p`i'resid, resid									//squared residuals for glm fit
	g p`i'residsq = (p`i'resid)^2
glm p`i'residsq b2.relig, family(gamma) link(log) 
predict p`i'residsqhat, mu 										// fitted variances, exp(Xb)
g LOGLIK=-.5*(ln(p`i'residsqhat)+(p`i'residsq/p`i'residsqhat))  // evaluating log likelihood
egen LL0 = sum(LOGLIK)											// summing log likelihood
display LL0
* Updating beta and lambda coefficients;
gen DLL=1														// initialize change in loglik;
while DLL > .00001 {
	drop p`i'resid
	reg p`i'_ i.site2 [aweight=1/p`i'residsqhat]		// WLS with variances as weights;
		drop p`i'residsqhat p`i'residsq LOGLIK
		predict p`i'resid, resid								// WLS resids;
		g p`i'residsq=(p`i'resid)^2								// squared resids for glm fit;
		estimates store p`i'BETA								// saving beta coefs;
	glm p`i'residsq b2.relig, family(gamma) link(log)			// gamma reg on log(r2);
		predict p`i'residsqhat, mu								// fitted variances, exp(Xb)
		estimates store p`i'LAMBDA								// saving lambda coefs;
	g LOGLIK=-.5*(ln(p`i'residsqhat)+(p`i'residsq/p`i'residsqhat))	// evaluating log likelihood;
	egen LLN = sum(LOGLIK)											// summing log likelihood;
		display LLN
	replace DLL=LLN-LL0											//  assess convergence;
	replace LL0=LLN
	drop LLN
}
drop p`i'resid p`i'residsqhat p`i'residsq LOGLIK LL0 DLL
}
esttab p21BETA p18BETA p19BETA p20BETA using vfr_firststage.tex, ///
		b(3) se(3) nostar replace label nobaselevels mtitles("Environment" "Gay Marriage" "Abortion" "Racism") 	///
		addnotes("Coefficients represent results from models regressing policy attitudes on a" ///
		"categorical variable for site. Evangelical Church 1 is the excluded category." ///
		"Standard errors in parentheses.") nonote // table of first stage regressions
esttab p21LAMBDA p18LAMBDA p19LAMBDA p20LAMBDA using vfr_secondstage.tex, ///
		b(3) se(3) nostar replace label nobaselevels mtitles("Environment" "Gay Marriage" "Abortion" "Racism") ///
		addnotes("Coefficients represent results from models regressing squared residuals from first" ///
		"stage on a categorical variable for type of site. Evangelical is the excluded" ///
		"category. Standard errors in parentheses.") nonote // table of second stage regressions

		